
#include <vector_types.h>
#include <math_functions.h>

struct Particula {
    float3 velocidade;
    float massaEspecifica;
    float pressao;
    float3 forcaResultanteChoque;
};

__device__ void atualizaPosicoesVBO(int indiceParticulaAtual, float3* posicoesParticulasT) {
    posicoesParticulasT[indiceParticulaAtual].x += 0.0001;
    posicoesParticulasT[indiceParticulaAtual].y += 0.0001;
    posicoesParticulasT[indiceParticulaAtual].z = posicoesParticulasT[indiceParticulaAtual].z;
}

__global__ void CalculosGPU(float3* posicoesParticulasT, float3* posicoesParticulasT2, int numeroParticulas, int resolucaoGridX, int resolucaoGridY, unsigned int* values, float deltaXGrid, float deltaYGrid, unsigned int* offSet, unsigned int* count, float h, Particula* particulasT, Particula* particulasT2, float massa, float constGas, float massaEspecificaInicial, float coeficienteViscosidade) {
    int thread = threadIdx.x;
    int block = blockIdx.x;
    int dim = blockDim.x;
    int indiceValue = dim * block + thread;
    int jAtual = 0;
    int iAtual = 0;
    int jMin = 0;
    int jMax = resolucaoGridX;
    int iMin = 0;
    int iMax = resolucaoGridY;
    int idCell;
    int kInicial;
    int kFinal;
    float3 forcaPressao = {0.0f, 0.0f, 0.0f};
    float3 forcaAceleracao = {0.0f, 0.0f, 0.0f};
    float3 forcaViscosidade = {0.0f, 0.0f, 0.0f};

    if (indiceValue < numeroParticulas) {
        int indiceParticulaAtual = values[indiceValue];
        jAtual = (posicoesParticulasT[indiceParticulaAtual].x / deltaXGrid);
        iAtual = (posicoesParticulasT[indiceParticulaAtual].y / deltaYGrid);

        //VERIFICAO DAS EXTREMIDADES.
        if (jAtual <= 0) {
            jMin = 0;
        } else {
            jMin = jAtual - 1;
        }
        if (jAtual >= (resolucaoGridX - 1)) {
            jMax = jAtual;
        } else {
            jMax = jAtual + 1;
        }
        if (iAtual <= 0) {
            iMin = 0;
        } else {
            iMin = iAtual - 1;
        }
        if (iAtual >= (resolucaoGridY - 1)) {
            iMax = (resolucaoGridY - 1);
        } else {
            iMax = iAtual + 1;
        }

        particulasT2[indiceParticulaAtual].massaEspecifica = 0;

        //CALCULO DA MASSA ESPECIFICA.
        /*for (int i = iMin; i <= iMax; i++) {
            for (int j = jMin; j <= jMax; j++) {
                float massaEsp = 0;
                idCell = (i * resolucaoGridX) + j;
                kInicial = offSet[idCell];
                kFinal = offSet[idCell] + count[idCell];
                for (int k = kInicial; k < kFinal; k++) {
                    int indiceParticulaVizinha = values[k];
                    if (indiceParticulaVizinha != indiceParticulaAtual) {
                        float distancia = sqrt((posicoesParticulasT[indiceParticulaAtual].x - posicoesParticulasT[indiceParticulaVizinha].x) * (posicoesParticulasT[indiceParticulaAtual].x - posicoesParticulasT[indiceParticulaVizinha].x) + (posicoesParticulasT[indiceParticulaAtual].y - posicoesParticulasT[indiceParticulaVizinha].y) * (posicoesParticulasT[indiceParticulaAtual].y - posicoesParticulasT[indiceParticulaVizinha].y));
                        if (distancia < h) {
                            //posicoesParticulasT[indiceParticulaAtual].x = (posicoesParticulasT2[indiceParticulaAtual].x + posicoesParticulasT2[indiceParticulaVizinha].x) * (distancia/h) /2;
                            //posicoesParticulasT[indiceParticulaAtual].y = (posicoesParticulasT2[indiceParticulaAtual].y + posicoesParticulasT2[indiceParticulaVizinha].y) * (distancia/h) /2;
                            //posicoesParticulasT[indiceParticulaAtual].z = posicoesParticulasT[indiceParticulaAtual].z;
                            massaEsp += massa * (315 / (64 * 3.14 * powf(h, 9)) * powf((h * h - distancia * distancia), 3));
                        }
                    }
                }
                particulasT2[indiceParticulaAtual].massaEspecifica += massaEsp;
            }
        }*/

        //CALCALA A PRESSAO NA PARTICULA.
        particulasT2[indiceParticulaAtual].pressao = constGas * (particulasT2[indiceParticulaAtual].massaEspecifica - massaEspecificaInicial);

        //CALCULA A FORÇA DE PRESSAO.
        for (int i = iMin; i <= iMax; i++) {
            for (int j = jMin; j <= jMax; j++) {
                float3 forcaPressaoTemp;
                forcaPressaoTemp.x = 0;
                forcaPressaoTemp.y = 0;
                forcaPressaoTemp.z = 0.0f;
                idCell = (i * resolucaoGridX) + j;
                float h2 = h*h;

                kInicial = offSet[idCell];
                kFinal = offSet[idCell] + count[idCell];

                for (int k = kInicial; k < kFinal; k++) {
                    int indiceParticulaVizinha = values[k];
                    if (indiceParticulaVizinha != indiceParticulaAtual) {
                        float distancia2 = (posicoesParticulasT[indiceParticulaAtual].x - posicoesParticulasT[indiceParticulaVizinha].x) * (posicoesParticulasT[indiceParticulaAtual].x - posicoesParticulasT[indiceParticulaVizinha].x) + (posicoesParticulasT[indiceParticulaAtual].y - posicoesParticulasT[indiceParticulaVizinha].y) * (posicoesParticulasT[indiceParticulaAtual].y - posicoesParticulasT[indiceParticulaVizinha].y);
                        if (distancia2 < h2) {
                            float distancia = sqrtf(distancia2);
                            float3 gradientePressao;
                            float mod;
                            mod = sqrtf(powf(posicoesParticulasT[indiceParticulaAtual].x - posicoesParticulasT[indiceParticulaVizinha].x, 2) +
                                    powf(posicoesParticulasT[indiceParticulaAtual].y - posicoesParticulasT[indiceParticulaVizinha].y, 2) +
                                    powf(posicoesParticulasT[indiceParticulaAtual].z - posicoesParticulasT[indiceParticulaVizinha].z, 2));
                            gradientePressao.x = (posicoesParticulasT[indiceParticulaVizinha].x - posicoesParticulasT[indiceParticulaAtual].x) / mod;
                            gradientePressao.y = (posicoesParticulasT[indiceParticulaVizinha].y - posicoesParticulasT[indiceParticulaAtual].y) / mod;
                            gradientePressao.z = (posicoesParticulasT[indiceParticulaVizinha].z - posicoesParticulasT[indiceParticulaAtual].z) / mod;
                            float kernelPress = (15 / (3.14 * powf(h, 6))) * (powf((h - distancia), 3));
                            float tempX = ((particulasT[indiceParticulaVizinha].pressao + particulasT[indiceParticulaAtual].pressao) / 2 * particulasT[indiceParticulaVizinha].massaEspecifica);
                            forcaPressaoTemp.x += massa * tempX * gradientePressao.x * kernelPress;
                            float tempY = ((particulasT[indiceParticulaVizinha].pressao + particulasT[indiceParticulaAtual].pressao) / 2 * particulasT[indiceParticulaVizinha].massaEspecifica);
                            forcaPressaoTemp.y += massa * tempY * gradientePressao.y * kernelPress;
                        }
                    }
                }
                forcaPressao.x += forcaPressaoTemp.x;
                forcaPressao.y += forcaPressaoTemp.y;
                forcaPressao.z += forcaPressaoTemp.z;
            }
        }
        forcaPressao.x *= -1;
        forcaPressao.y *= -1;
        forcaPressao.z *= -1;

        //CALCULA A FORÇA VISCOSIDADE.
        /*for (int i = iMin; i <= iMax; i++) {
            for (int j = jMin; j <= jMax; j++) {
                float3 forcaViscTemp;
                forcaViscTemp.x = 0;
                forcaViscTemp.y = 0;
                forcaViscTemp.z = 0;
                idCell = (i * resolucaoGridX) + j;

                kInicial = offSet[idCell];
                kFinal = offSet[idCell] + count[idCell];

                for (int k = kInicial; k < kFinal; k++) {
                    int indiceParticulaVizinha = values[k];
                    if (indiceParticulaVizinha != indiceParticulaAtual) {
                        float distancia = sqrtf((posicoesParticulasT[indiceParticulaAtual].x - posicoesParticulasT[indiceParticulaVizinha].x) * (posicoesParticulasT[indiceParticulaAtual].x - posicoesParticulasT[indiceParticulaVizinha].x) + (posicoesParticulasT[indiceParticulaAtual].y - posicoesParticulasT[indiceParticulaVizinha].y) * (posicoesParticulasT[indiceParticulaAtual].y - posicoesParticulasT[indiceParticulaVizinha].y));
                        if (distancia > 0 && distancia < h) {
                            float kernelVisc = (45 / (3.14 * powf(h, 6))) * (h - distancia);
                            forcaViscTemp.x += massa * ((particulasT[indiceParticulaVizinha].velocidade.x - particulasT[indiceParticulaAtual].velocidade.x) / 2 * particulasT[indiceParticulaVizinha].massaEspecifica) * kernelVisc;
                            forcaViscTemp.y += massa * ((particulasT[indiceParticulaVizinha].velocidade.y - particulasT[indiceParticulaAtual].velocidade.y) / 2 * particulasT[indiceParticulaVizinha].massaEspecifica) * kernelVisc;
                        }
                    }
                }
                forcaViscosidade.x += forcaViscTemp.x;
                forcaViscosidade.y += forcaViscTemp.y;
                forcaViscosidade.z += forcaViscTemp.z;
            }
        }*/
        forcaViscosidade.x *= coeficienteViscosidade;
        forcaViscosidade.y *= coeficienteViscosidade;
        forcaViscosidade.z *= coeficienteViscosidade;


        atualizaPosicoesVBO(indiceParticulaAtual, posicoesParticulasT);
    }
}

extern "C" void launch_CalculosGPU(unsigned int blocksPerGrid, unsigned int threadsPerBlock, float3* posicoesParticulasT, float3* posicoesParticulasT2, int numeroParticulas, int resolucaoGridX, int resolucaoGridY, unsigned int* values, float deltaXGrid, float deltaYGrid, unsigned int* offSet, unsigned int* count, float h, Particula* particulasT, Particula* particulasT2, float massa, float constGas, float massaEspecificaInicial, float coeficienteViscosidade) {
    // execute the kernel
    CalculosGPU << <blocksPerGrid, threadsPerBlock >> >(posicoesParticulasT, posicoesParticulasT2, numeroParticulas, resolucaoGridX, resolucaoGridY, values, deltaXGrid, deltaYGrid, offSet, count, h, particulasT, particulasT2, massa, constGas, massaEspecificaInicial, coeficienteViscosidade);
}
